function [DiagBlock, SampleList, RKTree, nDiagBlock] = ConstructHMatrix(N, nLevel)

solvetime = 0;
RKmulttime = 0;

DIRECT = 1;
INDIRECT = 0;

nWidthTop = N / 4;
% nWidthBottom = 32;
% nLevel = round( log2( nWidthTop / nWidthBottom ) + 1);

maxRank = zeros(nLevel, 1);
maxRank(1) = 20;
maxRank(2:end) = 10;

nSample = maxRank + 10;
% SVDCut = 1d-6;
SVDCut = zeros(nLevel, 1);
SVDCut(1) = 1e-6;
global svdcase;
switch svdcase
  case 1
    for i = 2 : nLevel
      SVDCut(i) = SVDCut(i-1) / 2;
    end
  case 2
    for i = 2 : nLevel
      SVDCut(i) = SVDCut(i-1) / sqrt(2);
    end
  case 3
    for i = 2 : nLevel
      SVDCut(i) = SVDCut(i-1);
    end
  case 4
    for i = 2 : nLevel
      SVDCut(i) = SVDCut(i-1) * sqrt(2);
    end
  case 5
    for i = 2 : nLevel
      SVDCut(i) = SVDCut(i-1) * 2;
    end
end 

SampleList = ConstructSampleList( N, nLevel );
RKTree = ConstructRKTree(N, nLevel, SampleList);

display('Constructing RKTree...');
if(1)
  for iLevel = 1 : nLevel
    display( iLevel )
    tic
    
    nWidth = RKTree{iLevel}.nWidth;
    nGroup = SampleList{iLevel}.nSource / ...
      SampleList{iLevel}.nParallel;
    TempR = randn( N, N, nSample(iLevel) );
    TempQ = cell(RKTree{iLevel}.nRKBlock, 1);

    % Preparing random matrices
    for i = 1 : RKTree{iLevel}.nRKBlock
      TempQ{i} = zeros( nWidth, nWidth, nSample(iLevel) );
    end

    % Applying the kernel to the random matrices
    for iGroup = 1 : nGroup

      TempMatApply = zeros( N, N, nSample(iLevel) );
      for iParallel = 1 : SampleList{iLevel}.nParallel
	iSource = (iGroup-1) * SampleList{iLevel}.nParallel + ...
	  iParallel;
	Source = SampleList{iLevel}.Source( :, iSource );
	TempMatApply( 1+(Source(1)-1)*nWidth : ...
	    Source(1)*nWidth, ...
	    1+(Source(2)-1)*nWidth : ...
	    Source(2)*nWidth, : ) = ...
	    TempR( 1+(Source(1)-1)*nWidth : ...
	    Source(1)*nWidth, ...
	    1+(Source(2)-1)*nWidth : ...
	    Source(2)*nWidth, : );
      end
      TempMatRes = KernelApply( N, reshape(...
	TempMatApply, N * N, nSample(iLevel) ) );
      TempMatRes = reshape( TempMatRes, N, N, nSample(iLevel) ); 
      TempMatRKRes = RKTreeApply( RKTree, TempMatApply, N, iLevel, nSample(iLevel) );
      TempMatRes = TempMatRes - TempMatRKRes;
        
      % Extracting the blocks for each target and perform SVD
      for iParallel = 1 : SampleList{iLevel}.nParallel
	iSource = (iGroup-1) * SampleList{iLevel}.nParallel + ...
	  iParallel;
	for iTarget = 1 : SampleList{iLevel}.nTarget
	  iRKBlock = SampleList{iLevel}.RKIndex(iTarget, iSource );
	  Target = SampleList{iLevel}.SourceTarget( :, ...
	    iTarget, iSource );
	  TempT = TempMatRes(...
	    1+(Target(1)-1)*nWidth : ...
	    Target(1)*nWidth,  ...
	    1+(Target(2)-1)*nWidth : ...
	    Target(2)*nWidth, : );
	  [TempU, TempS, TempW ] = svd( ...
	    reshape( TempT, nWidth * nWidth, nSample(iLevel) ), 'econ' );
	  nRankU = min( max(find( abs(diag(TempS)) ./ max(abs(diag(TempS))) > SVDCut(iLevel) )), ...
	    maxRank(iLevel) );
	  if( nRankU == maxRank(iLevel)  )
	    disp('max rank reached')
	  end
	  if( SampleList{iLevel}.SymFlag(iTarget, iSource ) == DIRECT )
	    RKTree{iLevel}.RKBlockU{ iRKBlock } = TempU(:, 1:nRankU );
	    TempQ{iRKBlock} = TempT;
	  else
	    RKTree{iLevel}.RKBlockV{ iRKBlock } = TempU(:, 1:nRankU );
	  end
	  clear TempT TempU TempS TempW;
	end % iTarget
      end % iParallel	

      clear TempMatApply TempMatRes TempMatRKRes;

    end % iGroup

    % Calculating Sigma and form low rank approximation
    for iRKBlock = 1 : RKTree{iLevel}.nRKBlock
      Source = RKTree{iLevel}.RKSource(:, :, iRKBlock);
      Target = RKTree{iLevel}.RKTarget(:, :, iRKBlock);
      TempR1 = reshape( ...
	TempR( Source(1,1):Source(1,2), Source(2,1):Source(2,2), : ), ...
	nWidth * nWidth, nSample(iLevel) );
      TempR2 = reshape( ...
	TempR( Target(1,1):Target(1,2), Target(2,1):Target(2,2), : ), ...
	nWidth * nWidth, nSample(iLevel) );

      TempSigma = pinv( transpose(TempR2) * ...
	RKTree{iLevel}.RKBlockU{iRKBlock} ) * ...
	transpose(TempR2) * reshape( ...
	TempQ{iRKBlock}, nWidth * nWidth, nSample(iLevel) ) * ...
	pinv( transpose(RKTree{iLevel}.RKBlockV{iRKBlock}) * ...
	TempR1 );
      sizeSigma = size(TempSigma);
      if( sizeSigma(1) > sizeSigma(2) )
	RKTree{iLevel}.RKBlockU{iRKBlock} = ...
	  RKTree{iLevel}.RKBlockU{iRKBlock} * ...
	  TempSigma;
	RKTree{iLevel}.RKBlockV{iRKBlock} = ...
	  RKTree{iLevel}.RKBlockV{iRKBlock};
      else
	RKTree{iLevel}.RKBlockU{iRKBlock} = ...
	  RKTree{iLevel}.RKBlockU{iRKBlock};
	RKTree{iLevel}.RKBlockV{iRKBlock} = ...
	  RKTree{iLevel}.RKBlockV{iRKBlock} * ...
	  transpose(TempSigma);
      end
      clear TempR1 TempR2 TempSigma;
    end % iRKBlock
    
    clear TempR TempQ;

    toc
  end % iLevel
end

% Extract diag part

display('Extracting diagonal...');
tic

% nWidth is doubled because we cannot sample two adjacent blocks.
% nBlockPre = round( N / nWidth );
nWidthPre = nWidth;
nWidth = nWidth * 4;
nBlockCur = round( N / nWidth );

if( 0 ) % when there is access to large memory
  nDiagBlock = 4*4*nBlockCur*nBlockCur;
  DiagBlock = cell(nDiagBlock, 1);
  for i = 1 : nDiagBlock
    DiagBlock{i} = struct('indSource', [], ...
      'indTarget', [], ...
      'Mat', [] );
  end
  nSampleDiag = nWidth * nWidth;

  disp('Generating pattern');
  TempMatApply = zeros( N, N, nSampleDiag );
  for iSourceX = 1 : nWidth
  for iSourceY = 1 : nWidth
      SourceYInd = iSourceY : nWidth : N;
      SourceXInd = iSourceX : nWidth : N;
      ind = iSourceY + (iSourceX-1) * nWidth;
      TempMatApply(SourceYInd, SourceXInd, ind) = 1;
    end
  end

  disp('Calculating');
  TempMatRes = KernelApply( N, reshape(...
    TempMatApply, N * N, nSampleDiag) );
  TempMatRes = reshape( TempMatRes, N, N, nSampleDiag ); 
  TempMatRKRes = RKTreeApply( RKTree, TempMatApply, N, nLevel+1, nSampleDiag );
  TempMatRes = TempMatRes - TempMatRKRes;
  TempMatRes = reshape( TempMatRes, N*N, nSampleDiag );

  disp('Redistributing');
  cDiagBlock = 0;
  for j = 1 : 4
    for i = 1 : 4
      iTempSourceArr = ( ((i-1)*nWidthPre + 1) : i*nWidthPre )';
      jTempSourceArr = ( ((j-1)*nWidthPre + 1) : j*nWidthPre );
      indTempSource = repmat( (jTempSourceArr-1)*nWidth, length(iTempSourceArr), 1 ) + ...
	repmat( iTempSourceArr, 1, length(jTempSourceArr) );
      indTempSource = reshape(indTempSource, [], 1);
      for k = 1 : nBlockCur
	for l = 1 : nBlockCur
	  iSource = (k-1)*nWidth + (i-1)*nWidthPre + 1;
	  jSource = (l-1)*nWidth + (j-1)*nWidthPre + 1;
	  iSourceArr = (iSource:iSource+nWidthPre-1)';
	  jSourceArr = (jSource:jSource+nWidthPre-1);
	  indSource = repmat( (jSourceArr-1)*N, length(iSourceArr), 1 ) + ...
	    repmat( iSourceArr, 1, length(jSourceArr) );
	  indSource = reshape(indSource, [], 1);
	  iSourceBlk = (k-1)*4 + i;
	  jSourceBlk = (l-1)*4 + j;
	  iTarget1 = (iSourceBlk-2) * nWidthPre + 1;
	  iTarget2 = (iSourceBlk+1) * nWidthPre;
	  jTarget1 = (jSourceBlk-2) * nWidthPre + 1;
	  jTarget2 = (jSourceBlk+1) * nWidthPre;
	  iTargetArr = (iTarget1:iTarget2)';
	  jTargetArr = jTarget1:jTarget2;
	  iTargetArr = mod( iTargetArr - 1, N ) + 1;
	  jTargetArr = mod( jTargetArr - 1, N ) + 1;
	  indTarget = repmat( (jTargetArr-1)*N, length(iTargetArr), 1 ) + ...
	    repmat( iTargetArr, 1, length(jTargetArr) );
	  indTarget = reshape(indTarget, [], 1);
	  cDiagBlock = cDiagBlock + 1;
	  DiagBlock{cDiagBlock}.indSource = indSource;
	  DiagBlock{cDiagBlock}.indTarget = indTarget;
	  DiagBlock{cDiagBlock}.Mat = ...
	    (TempMatRes(indTarget, indTempSource)).';
	end
      end
    end
  end

end

if( 1 ) % when memory is not enough
  % The first dimension: i, k, or Y
  % The second dimension: j, l, or X
  %
  nDiagBlock = 4*4*nBlockCur*nBlockCur;
  DiagBlock = cell(nDiagBlock, 1);
  for i = 1 : nDiagBlock
    DiagBlock{i} = struct('indSource', [], ...
      'indTarget', [], ...
      'Mat', [] );
  end
  nSampleDiag = nWidth * nWidth;

  cDiagBlock = 0;
  for j = 1 : 4
    for i = 1 : 4
      for k = 1 : nBlockCur
	for l = 1 : nBlockCur
	  iSource = (k-1)*nWidth + (i-1)*nWidthPre + 1;
	  jSource = (l-1)*nWidth + (j-1)*nWidthPre + 1;
	  iSourceArr = (iSource:iSource+nWidthPre-1)';
	  jSourceArr = (jSource:jSource+nWidthPre-1);
	  indSource = repmat( (jSourceArr-1)*N, length(iSourceArr), 1 ) + ...
	    repmat( iSourceArr, 1, length(jSourceArr) );
	  indSource = reshape(indSource, [], 1);
	  iSourceBlk = (k-1)*4 + i;
	  jSourceBlk = (l-1)*4 + j;
	  iTarget1 = (iSourceBlk-2) * nWidthPre + 1;
	  iTarget2 = (iSourceBlk+1) * nWidthPre;
	  jTarget1 = (jSourceBlk-2) * nWidthPre + 1;
	  jTarget2 = (jSourceBlk+1) * nWidthPre;
	  iTargetArr = (iTarget1:iTarget2)';
	  jTargetArr = jTarget1:jTarget2;
	  iTargetArr = mod( iTargetArr - 1, N ) + 1;
	  jTargetArr = mod( jTargetArr - 1, N ) + 1;
	  indTarget = repmat( (jTargetArr-1)*N, length(iTargetArr), 1 ) + ...
	    repmat( iTargetArr, 1, length(jTargetArr) );
	  indTarget = reshape(indTarget, [], 1);
	  cDiagBlock = cDiagBlock + 1;
	  DiagBlock{cDiagBlock}.indSource = indSource;
	  DiagBlock{cDiagBlock}.indTarget = indTarget;
	end
      end
    end
  end

  nVecMax = 400;

  TempMatApply = zeros( N, N, nVecMax );
  TempMatRes = zeros(N, N, nVecMax);
  TempMatRKRes = zeros(N, N, nVecMax);
  cVec = 0;
  arrSourceY = zeros(nVecMax);
  arrSourceX = zeros(nVecMax);
  for iSourceX = 1 : nWidth
    for iSourceY = 1 : nWidth
     
      cVec = cVec + 1;
      arrSourceY(cVec) = iSourceY;
      arrSourceX(cVec) = iSourceX;
      SourceYInd = iSourceY : nWidth : N;
      SourceXInd = iSourceX : nWidth : N;
      TempMatApply(SourceYInd, SourceXInd, cVec) = 1;
      if( (cVec < nVecMax) && ((iSourceX<nWidth) || (iSourceY<nWidth)) )
	continue;
      end
     
      TempMatRes = KernelApply( N, reshape(...
	TempMatApply, N * N, nVecMax) );
      TempMatRes = reshape( TempMatRes, N, N, nVecMax ); 
      TempMatRKRes = RKTreeApply( RKTree, TempMatApply, N, nLevel+1, nVecMax);
      TempMatRes = TempMatRes - TempMatRKRes;
      TempMatRes = reshape( TempMatRes, N*N, nVecMax);

      for iVec = 1 : cVec
	YSourceBlk = floor((arrSourceY(iVec)-1)/nWidthPre) + 1;
	XSourceBlk = floor((arrSourceX(iVec)-1)/nWidthPre) + 1;
	indSourceBlk = (XSourceBlk-1)*4 + YSourceBlk;
	YBlock = mod( arrSourceY(iVec)-1, nWidthPre ) + 1;
	XBlock = mod( arrSourceX(iVec)-1, nWidthPre ) + 1;
	indBlock = (XBlock-1)*nWidthPre + YBlock;
	for l = 1 : nBlockCur
	  for k = 1 : nBlockCur
	    cDiagBlock = nBlockCur*nBlockCur*(indSourceBlk-1) + ...
	      (k-1)*nBlockCur + l;
	    DiagBlock{cDiagBlock}.Mat(indBlock, :) = ...
	      (TempMatRes(DiagBlock{cDiagBlock}.indTarget, iVec)).';
	  end
	end
      end
	
      cVec = 0;
      TempMatApply = zeros( N, N, nVecMax );
      TempMatRes = zeros(N, N, nVecMax);
      TempMatRKRes = zeros(N, N, nVecMax);

    end
  end


end

toc
whos
